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Abstract 

We propose numerical simulations of the Ashkin- Teller model as a foil for 
theoretical techniques for studying very weakly first-order phase transitions in 
three dimensions. The Ashkin- Teller model is a simple two-spin model whose 
parameters can be adjusted so that it has an arbitrarily weakly first-order 
phase transition. In this limit, there are quantities characterizing the first- 
order transition which are universal: we measure the relative discontinuity 
of the specific heat, the correlation length, and the susceptibility across the 
transition by Monte Carlo simulation. 
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I. INTRODUCTION 



One of the modem scenarios for explaining the baryon asymmetry of the universe de- 
pends on the nature of the electroweak phase transition in the early universe. The scenario 
requires the transition to be first-order, but in some cases of interest the transition is suffi- 
ciently weak and near-critical that simple perturbative expansions around mean field theory 
are not well-behaved. The study of this cosmological application has spawned a general 
renewal of interest in techniques for studying critical or near-critical phase transitions: tech- 
niques that have a two decade history in condensed matter physics. The Ising model is the 
canonical example of a three-dimensional system with a second-order phase transition and 
has been well studied both numerically and theoretically. There has not, however, been a 
comparable amount of attention paid to any canonical example of weak, near-critical first- 
order transitions. A simple two-spin generalization of the Ising model that can provide such 
a canonical example — a testbed for theoretical techniques for treating very weakly first-order 
transitions in three dimensions — is the Ashkin- Teller model. By tuning parameters in this 
model, one can obtain first-order transitions that are arbitrarily weak. 

A full discussion of our motivation for studying this model, and its similarities and 
dissimilarities with the electroweak phase transition, is given in ref. In the present 
paper, we will study the relative discontinuity of various physical quantities across the first- 
order transition. In particular, we measure the ratios C+/C_, ^+/^_, and X+/X- where C+, 

x+ ci-re the specific heat, correlation length, and susceptibility in the disordered (high- 
temperature) phase and C_, and ^'^^ the same in the ordered (low-temperature) 
phase. These ratios are universal in the limit that the first-order transition is arbitrarily 
weak and provide a set of tests against which to measure theoretical techniques for studying 
weakly first-order transitions. In particular, a comparison of our numerical results against 
the predictions of e-expansion methods may be found in ref. 

In the remainder of this introduction, we review the Ashkin- Teller model and then present 
our final results. We will review the parameter x of the Ashkin- Teller model whose *-0"'" 
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limit yields arbitrarily weakly first-order transitions. In sec. |T|, we give a broad overview of 
our method for making measurements in the two phases and for determining the transition 
temperature. In sec. |IT|, we present our measured ratios for an assortment of values of x] 
discuss how C±, and x± should scale at small x; explain how finite-x corrections to our 
X— s>0+ ratios should scale with x; and discuss our procedure for extracting the x— S'O"'" limit. 
In sec, we discuss the details that went into the individual measurements at each x such 
as our procedure for measuring susceptibilities and correlation lengths, our assessment of 
finite volume errors, and uncertainties in the transition temperature and their effect on our 
measurements. 



A. Review of the Ashkin- Teller Model 

The (symmetric) Ashkin- Teller model ^ is a system with two Ising spins Sj and ti per 
lattice site z, with the nearest-neighbor interaction 

[5H = -f5^{siSj + titj + xsiUsjtj) , Si,tj = ±l. (1.1) 

Special cases of interest include x=0, which corresponds to two decoupled Ising models, and 
x=l, which can be rewritten as 



fiH{x=l) = -A(3j2{Ss.sA.t, - i) (1.2) 

and is equivalent to the 4-state Potts model [Q. The phase diagram of the model in three 
dimensions on a simple cubic lattice is sketched in fig. |I| 0. The portion of this diagram 
we will focus on is the neighborhood of the Ising tri-critical point at x=0. By studying 
arbitrarily small but positive values of x, we can study arbitrarily weak first-order phase 
transitions. Our goal will be to extract ratios such as 

lim — = lim , (1.3) 

x^o+ (7_ x~^o+ hmp^^+ C 

where Pt = Pt{x) is the inverse transition temperature for a given value of x. 
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FIG. 1. Believed phase diagram of three-dimensional Ashkin- Teller model on a simple cubic lattice, 
taken from ref. where it was extracted from series analysis and Monte Carlo data. Dashed and solid 
lines indicate first and second-order transitions respectively. Dotted lines indicate cases where the nature of 
the transition has not been unambiguously determined. The phases are labeled disordered ((s) — {t) — 0); 
Baxter (ferromagnetic with (s), (i), and (st) non-zero); "(si)" (where {st) is ferromagnetically ordered but 
(s) = (t) = 0); "(si)AF" (the same but anti-ferromagnetically ordered); and "(s)" (where either (s) or {t) is 
ferromagnetically ordered but the other is not and (st) — 0.) 
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For generic values of the parameter x, the internal symmetry group of the model is that 

of a square, D4, acting on the spin states {s,t) = (++), (H — ), ( ), ( — h). The line x=0 is 

a line of enhanced symmetry: the translation symmetry decouples for the s and t spins.Q (A 
simple continuum model with the same long-distance degrees of freedom and symmetries is 
the cubic anisotropy model, a two-scalar field theory discussed in refs. 



B. Final results 

Our final results for the universal ratios are[| 

0.071(8), (1.4) 

1.6(2), (1.5) 
4.0(6), (1.6) 
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where the error estimate on the last result should be taken with a grain of salt (see sec. [Ill C 
for details). The susceptibility x is formally defined in the infinite- volume limit by adding 
a source term — Xli h ■ Sj to PH, where S = (s, t), and then taking 

X = ^limE-^(5.). (1.7) 

2 h^O ^ dha 

As an incidental consequence of our work, we also have what is, to our knowledge, the 
best measurement of the transition temperature of the 4-state Potts model (a;=l) in three 
dimensions: 

0.157154(4), simple cubic; 
/3t(4-state Potts) = { (1.8) 

0.113752(5), body-centered cubic; 



^ This makes the Ising point a particularly attractive muhi-critical point to study because we know it's 
value of X in advance. Other multi-critical points would first require a numerical search for the corresponding 
multi-critical value Xc of x before we could proceed to numerically extract the limit x—^Xc of ratios analogous 



to (1.3) 



^ Due to an improvement in our method of error analysis prior to publication, the numbers given here are 
slightly, but not significantly, different than those originally reported in ref. . 
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where j3 is normalized according to ( p,.2| ). For comparison, earlier values determined by 
Monte Carlo and series analysis are given in Table |. The series estimate of ref. |Tl| for 
body-centered cubic (BCC) lattices is off by roughly thrice their estimated error. 
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simple cubic 
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Monte Carlo, 


ref. i 


0.158 




Monte Carlo, ref. lU 


~ 0.15694 




series, ref. 
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0.162(2) 


0.1176(13) 


series, ref. 


12 




0.161(3) 





TABLE I. Summary of previous Monte Carlo and series analysis results for j3t in the 4-state Potts 
model (a^=l) in three dimensions. The Monte Carlo references do not quote error estimates. Note that our 
convention (^]^) for (3 differs from the usual one by a factor of 4. 



II. OUTLINE OF CALCULATIONAL METHOD 

Our lattices are simple cubic with helical boundary conditions (periodic except for a 
twist) and range in sizes up to 480 x 120^. In principle, it would be nice to verify the 
universality of a;^0^ ratios such as C+/C_ by repeating the calculation on another lattice 
type, such as body-centered cubic. However, we have only made limited exploration of BCC 
lattices, and our data for this case is relegated to Appendix 0. 

We interlace two different update algorithms: a simple heatbath algorithm and a cluster- 
update algorithm.0 The cluster-update algorithm is a simple generalization of the usual 
cluster-update algorithms for 0{n) systems [|l^] and is described in Appendix 0. 

Our basic method for extracting our results is best outlined with a specific example: 
the specific heat ratio C+/C_. The first step is to measure C+/C_ at fixed values of x, so 



But some of our data was generated by heatbath only or by cluster only. 
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FIG. 2. Specific heat vs. (3 at fixed x for (a) infinite volume, (b) finite volume, (c) finite-time numerical 
simulations in large volume. The (5-function in (a) at /3t represents the latent heat. In all cases, [3t represents 
the true, infinite-volume transition temperature. 

that we can later extract the s^O^ hmit. Fig. ^ quaUtatively sketches the dependence 
of specific heat on inverse temperature (3 for fixed x. In finite volume, the latent heat 5- 
function at the transition temperature broadens out into a finite-width peak as in fig. ^b. 
If we actually measured this peak in numerical simulations, we would have to attack the 
problem of how to extract Cj^ and C_ from beneath the peak. However, the lattices we 
work on are large enough that the mixing time between the two phases at the transition is 
very long compared to the duration of our simulations.^ Our results therefore have the form 
of fig. there is a range of temperatures, around the transition temperature, in which 
our simulation sits in either the ordered or disordered phase depending on whether we start 
with ordered or disordered initial conditions. The problem of extracting C+ and C_ now 
becomes the problem of determining the inverse transition temperature (3^. 

We determine the transition temperature by working on asymmetric lattices of size L x 
T X T with L > T. We divide the L dimension in half, and we start the lattice in the 



This is not an accident. The cleanest way to eliminate systematic errors from finite- volume effects is to 
work in volumes V large compared to the correlation volume V^. But, at the transition temperature, the 
mixing time grows exponentially with V/V^. 
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FIG. 3. Initial conditions used on asymmetric lattices for determining the transition temperature. 

ordered phase in one half and the disordered phase in the other, as depicted in fig. |^. Next 
we evolve the system with our Monte Carlo update procedure. The domain walls will feel a 
net pressure to expand the phase with the lowest free energy. By tracking whether the system 
eventually evolves into the ordered or disordered phase, and finding the value of (3 where the 
favored phase changes, one determines the transition temperature. The story is a little more 
complicated, however. Exactly at the transition temperature, neither phase is favored and 
the domain walls will random walk until, randomly, they collapse the system into one phase 
or the other. Fig. ^ shows an example from our simulations of collapsing into either phase at 
the transition temperature. Slightly away from the transition temperature, there will be a 
slight bias to the random walks, but there will still be some chance of ending in the disfavored 
phase. The probability of ending up in one particular phase, say the ordered one, therefore 
has a dependence on temperature like that sketched in fig. It can be modeled as the 
solution to a classic problem from probability theory — the Gambler's Ruin problem — which 
is described in Appendix p. Our technique for determining (3^ is therefore to make multiple, 
independent runs at each temperature, in order to numerically extract the curve of fig. |^a, 
and then to fit for /3t. There is one more complication: even if the transverse dimension 
T is large compared to the correlation length, there are systematic errors in this procedure 
that have only power-law fall-off with increasing longitudinal dimension L, as depicted in 
fig. To get the best estimate of the transition temperature, we make the additional step 
of numerically extracting the L—*oo limit. A model for the finite L corrections is discussed 
in Appendix |B] and sec. |IVD . 



All of our calculations were carried out on a handful of SGI Indy R4600 workstations. 



9 



-2 



a) 

i3X) -3 

d 



-4 



disordered 

i/' 



ordered 



500 1000 

Monte Carlo time 



1500 



FIG. 4. Examples of evolution of energy density with Monte Carlo time starting from the mixed phase 
of fig. ^. The data is shown for on a 160 x 40^ lattice near the transition temperature {/3 = 0.180263). 




FIG. 5. (a) Probability of ending in the ordered phase vs. (3 when starting with mixed initial conditions; 
(b) the same but including finite L effects. The narrower curves correspond to larger L. 
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taking roughly 2 to 45 mins. per energy decorrelation time (depending on x) on an 80^ 
lattice at the transition. 



III. ANALYSIS OF DATA AND x^O 



Table |I| and fig. ^ summarize our results for various ratios as a function of x on simple 
cubic lattices; table |T| shows in more detail the individual quantities that went into deter- 
mining the ratios. We will explain in more detail our determination of these quantities in 
sec. 0. In this section, we will focus on understanding the small x behavior of the system 
and the extraction of the x —>■ limits of the ratios. The errors quoted for our data are 
mostly statistical but include the systematic error from our uncertainty in the transition 
temperatures (3f We estimate finite size effects to be no larger than our quoted total errors 
except possibly for our lowest value x=0.3 of x, where we do not have a very good estimate. 
For finite correlation length, the correlation length measured along edges of the lattice need 
not be exactly the same as that measured along diagonals, and so we have listed these cases 
separately. The x—>-0^ ratio ^+/^- in the two cases, however, should be the same. 




1.0 0.1228(19) 1.108(16)^ 

0.8 0.122(3) 1.191(24)^ 

0.6 0.109(3) 1.307(24)^ 

0.5 0.103(3) 1.394(21)^ 

0.3 0.088(8)^ 1.45(7)^'*^ 



1.062(12)^ 2.80(6) 
1.235(16)^ 3.77(9) 
1.333(20)^ 4.20(13) 



1.342(17)^ 
1.41(8)^''^ 



4.28(15) 
3.8(5)^ 



TABLE II. Summary of relative discontinuities of specific heat, correlation length (along lattice edges 
and diagonals), and susceptibility vs. x. Errors include statistical errors and systematic errors due to the 
uncertainty in the determination of transition temperatures. Finite volume errors are not included and are 
estimated to be no larger than the total errors quoted above, except possibly for (a) the case x=0.3 where 
we do not have a good estimate of them. (See sec. |^ for details.) (b) Errors for some correlation lengths 
are quite likely underestimated (see sec. |IVB for details) . 
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FIG. 7. The correlation length ^'^^sc (measured along lattice edges) at the transition vs. x for the 
disordered phase (squares) and ordered phase (diamonds). 

Our lowest value x = 0.3 of x is small in the sense that the transition is relatively weak: 
the correlation length at the transition is order 20 (see table HI). Indeed, it is precisely 
the rapid growth of correlation length with decreasing x shown in fig. |^ (combined with the 
requirement that lattices be large compared to ^ to avoid transitions between the phases) 
that has prevented us from simulating even smaller x. However, despite the fact that x = 0.3 
is "small," the quality of our data in fig. |^ is such that, when attempting to extract the x— >0 
limits, we would clearly benefit tremendously from knowing a priori how the correction to 
the x — limit should scale for small x. 

We shall first discuss how the individual quantities C±, and x± scale with x. Then 
we turn to the scaling of the corrections to x^O limit of ratios. Finally, armed with this 
analysis, we shall fit the data of fig. || as best we can. 
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1 


0.157154(4) 


-2.2734(3) 


1.663(13) 


2.878(24) 


2.985(14) 


33.1(3) 






-4.6009(23) 


13.54(18) 


2.60(3) 


2.81(3) 


11.84(22) 


0.8 


0.168149(3) 


-2.1552(4) 


1.836(24) 


3.62(3) 


3.74(3) 


51.9(9) 






-4.1726(23) 


15.1(4) 


3.04(6) 


3.03(3) 


13.78(20) 


0.6 


0.180272(3) 


-2.0942(4) 


2.33(3) 


5.32(4) 


5.47(5) 


111(3) 






-3.5359(23) 


21.3(4) 


4.07(7) 


4.10(5) 


26.3(3) 


0.5 


0.186750(4) 


-2.0797(3) 


2.72(3) 


7.28(5) 


7.29(4) 


203(3) 






-3.145(3) 


26.4(6) 


5.23(7) 


5.43(6) 


47.3(1.3) 


0.3 


0.2003659(15) 


-2.0658(4) 


5.16(21) 


22.1(7) 


21.6(8) 


1730(180) 






-2.362(3) 


58(4) 


15.2(5) 


15.3(6) 


460(30) 





0.221652(3)^^ 













TABLE III. Summary of inverse transition temperature /3t, energy density e±, specific heat density C±, 
correlation length and susceptibility x± a function of a; on a simple cubic lattice. + and — denote 
the disordered and ordered phases, respectively, at the transition temperature. For some x, the errors on 



the correlation lengths are likely underestimated (see LVB for details). The x~0 Ising model (3t (a) is taken 
from rcf. iH . 



A. The Problem with Crossover Exponents 

As we shall review below, one already knows how various dimensionful quantities, such 
as the correlation lengths should diverge as x^O. The divergence is characterized by 
crossover (or "tricritical" ) exponents, e.g. 

^± ~ x-\ x^O. (3.1) 

In the case at hand, crossover exponents such as y may be determined from knowledge of 
Ising model critical exponents. In the Ising model, the scaling dimension of the nearest- 
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neighbor interaction in the Hamiltonian density is 

Ida 

dim[siSi+e\ = "^~=2^2z/' ^ 
where d=3 is the dimension of space. Now consider the scahng dimension of the interaction 
{st)i{st)ij^e in the Ashkin- Teller model. For the purpose of understanding the crossover 
behavior at very small x, we can ignore the effect of the interaction itself on its scaling and 
instead consider the scaling dimension of the operator in the x=0 limit. Then the operator 
factorizes: 

D = dim[(st)i(st)j+e] = dim[siSj+e] + dim[titi+e] = ^~ ~- (3-3) 

Deviation from Ising behavior will occur once this operator becomes important, which in 
our case means that the scale ^ of the first-order transition is related to the interaction's 
coefficient x by 

X ~ ^^-'^ ~ ^-"/^ . (3.4) 

Using Ising scaling laws to get the power-law relationship of with other quantities, one 
gets 



C ^—5.7(3) 


(3.5a) 




(3.5b) 




(3.5c) 


model exponents 15 




z/ ~ 0.6300(15) , 7 - 1.2405(15). 


(3.6) 



This result means, in the small x limit, that the correlation length and susceptibility 
should grow by factors of roughly 40-60 and 1500-3500 respectively when x is reduced by a 
factor of two! There is clearly no sign of such strong x dependence in the data of table ^ 
and fig. 1^. A natural question now arises: Does this discrepancy indicate failure to reach the 
small X region of the Ashkin- Teller model and so make our x — limits for ratios suspect? 
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To address this question, one must first note that the origin of the strong x dependence 
in ( p.5| a-b) is the small value of the specific-heat exponent a in the Ising model. It is useful 
to formally consider what would have happened if a were arbitrarily small, or even zero. In 
Appendix 0, we consider a simple model for the renormalization group fiow where one can 
treat this explicitly. In the case a=0, the relationship ( |3.4| ) generically becomes a: ~ 1/ In(^), 
and one gets 

(3.7a) 
(3.7b) 
(3.7c) 

1 , (3.7d) 

where k is some constant. For arbitrarily small but non-zero a, one finds exponential scaling 
( p.7| a-b) of ^ and x for moderately small x, but the power- law scaling form ( |3.5| a-b) takes 
over when x drops below 

ak 

This threshold has no effect, however, on the extraction of the limits of ratios such as ^+/^_. 
(See Appendix for details.) So, though the formal limit x <^ 0{a) is required to see correct 
scaling of ^ and Xi the weaker limit x <^ 0(1) is adequate for extracting the dimensionless 
ratios of interest. Our failure to see the correct scaling of ^ and x our numerical data 
should not, therefore, cause concern. Roughly, we expect that x should be small enough for 
computing ratios if the correlation length is large compared to the lattice spacing. 

We can test the above assertions by examining relationships that don't depend on a and 
which should be satisfied for both the different scaling behaviors ( |3.5| ) and ( p.7| ). One such 
relationship is that S,± ~ Xi±'^'i therefore 

lim = - . 3.9 

x^o+ In x± 7 



Fig. ^ shows this ratio of logarithms for our data and the theoretical limit ( ^.9] ). The data 
appears consistent with the limit. Another test is to check whether xC± approaches a 
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FIG. 8. A test of scaling: ln(^®'^se^^ jj^^'-^'j yg_ ^ fgj. ^jjg ordered phase (diamonds) and disordered 
phase (circles). The cross marks the theoretical value v/j for the a;^0 limit. The corresponding graph for 
\n{(^'^^^^) / \n{x) is visually almost identical. 

constant as Fig. ^ shows xC+ for our data, which looks reasonably good. Fig. ^ 

shows xC-, which is reasonable except that the lowest x point is a bit high. (We shall see 
below that the approach to the x—>-0 limit should be linear at small x.) 



B. Corrections to Scaling 

In order to extrapolate x— >0 limits for our ratios in fig. ^, it helps tremendously to know 
how the corrections to those limits scale with x. For simplicity, let us first ignore the small a 
issue. We have discussed in the previous section how, for very small x, the system has Ising 
behavior for distance scales up to order ^ ~ x~'^^°'. Corrections to scaling arise because there 
are irrelevant operators which have not quite scaled away for finite (but large) ^. However, 
since the scaling of operators as ^ ^ cxd is determined by Ising behavior, we can extract the 
dimension of the most important such operator from well-known results in the Ising model: 
relative corrections to scaling behavior scale as where u is the Ising model exponent 

m 
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FIG. 9. A test of scaling: (a) xC+ and (b) xC- vs. x. 



Translating ^ to x using 



uj ~ 0.79(3) 



we then have, for example, 



X± 
X- 



0(1) + 0(a;' 



LL)v/a\ 



Similarly, 



^ ~0(l) + 0(x' 



(3.10) 



(3.11) 
(3.12) 



(3.13) 



The specific heat ratio is slightly different. In addition to the divergent 0{^°'/'^) = 0{x~^) 
contribution from long-distance modes, the specific heat receives a direct contribution from 
short- distance modes that is analytic in ^ and is therefore 0{^^) = 0{x^). This is much 
more important than the relative contribution discussed above and gives 



C± ~ 



o(r/^) + o(e°) 

0(l) + 0(x). 



(3.14) 
(3.15) 



So far, we have discussed the correction to scahng in the hmit of arbitrarily small x 
where everything scales with the correct Ising model exponents. But, as discussed in the 
previous section, the values of x we actually simulate are not small enough to reproduce 
relationships involving the small Ising exponent a. The scaling ( |3.15| ) of corrections to 



C+ / C_ is nonetheless in good shape because it does not depend on a and would hold even 
if a were zero. The scaling of corrections to other ratios can also be cast in a form that does 
not depend on a: 

^^0{l) + 0{x7"'''). (3.16a) 

X- 

|^~0(1) + 0(CT. (3.16b) 

( 3.15| ) and (|3.16|) are the forms we fit our data to in order to extract the x— >0 limits of the 
ratios. 

C. Extraction of x— >0 limits from data 

Now that we know how the corrections are supposed to scale, we can attempt to fit our 
data to the appropriate form. C+/C_ is already plotted against the correct variable x in 
fig. ^ for a linear fit. In fig. 0, we show x+IX- and i+H- plotted against x+"'^'^ a^iid ^'^'^ 
respectively, so that the fit should again be linear. Our procedure is to fit the largest set of 
points, working from the lowest x up, that yields a reasonable chi-squared. 

For C+/C_, fitting the points x < 0.8 yields a freakishly high confidence level of 98% and 
produces our final result, eq. (|1.4| ), for the x— >0 limit. Adding the point x = 1.0 decreases 
the confidence level to 18% and would change the limit to 0.084(4). Because of its large 
uncertainty, our x=0.3 value has an almost negligible effect on the fit and our final result. 

Next consider x+/X-- Fitting the three points x < 0.6 yields the best fit line (34% 
confidence level) shown in fig. |10|a and our final result, eq. (|1.6|) , for the x— >0 limit. The 



systematic uncertainty in the values of the Ising critical exponents is included in our error 
estimate. Attempting to add the x = 0.8 point gives lim(x_|_/x_) = 4.8(3) with a fairly small 
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FIG. 10. Ratios vs. an appropriate scaling variable for corrections to the x~^0 limit: (a) x+/x- "^s. 
and (b,c) vs. ^^'^ for edges and diagonals. Statistical uncertainties in the ordinates of 

the data points are too small to be seen. We have not shown for each data point the larger systematic 
uncertainties due to uncertainties in the Ising exponents (dominated by lo). The solid line is our best fit, 
and the x—0 value is our extrapolation. The dotted lines show how our best fit changes as uj is varied over 
the uncertainty indicated in ( ^.10 ), but this change is negligible in (a). 
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FIG. 11. ^ ratios vs. instead of ^ 
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confidence level of 11%. We have parameterized corrections ( p.l6|) to scaling by x^^^"^ ■ An 
alternative choice would be x^^^^ ^ and we have used x+^^^' simply because it is more 
straightforward to measure. If the fit is made against x^^^^ instead of 'x^^^'^ ^ the result 
for the limiting ratio is 4.0(5) and is consistent with the previous result. In any case, the 
situation is somewhat unsatisfactory to the extent that (1) we have fit only three points, 
and (2) one of those points is our lowest x value, x = 0.3, where we are less confident about 
finite volume errors (see sec. p.V C] ). 



In sec. |IV B| , we explain that our errors on extracting the correlation lengths may be 
underestimated. Our procedure for estimating our final error will simply be to interpolate 
the X— s>0 ratio in different ways that should in principle be equivalent. For ^'^^^ / ^'^^'^ , we can 
obtain a reasonable (56% confidence level) fit vs. (^^^so^-oj ^j^]-^ qH q^j- (iata points, shown 
in fig. p!0|b; for the diagonal ratio, we need to drop the x=1.0 point from the fit (improving 
the confidence level from O(10~^%) to 36%). The results are 

t-edge fdiag 

lim = 1.64(3) , lim = 1.51(4) . (3.17a) 



Alternate interpolations using are shown in fig. ^ and give 

t-edge fdiag 

lim = 1.72(5) , lim ^ = 1.53(5) . (3.17b) 

The wide discrepancy of values in (|3.17aD and (|3.17b| ) should be taken as a reflection of our 
systematic errors in determining the correlation lengths and in extracting the >0 limit. 
Our final result (|1.5| ) has been chosen to span all of the above extrapolations. 

Clearly it would be useful to have more good-quality data at small x, and, in particular, 
X = 0.4 suggests itself as a good candidate for future simulation. We estimate that x = 0.4 
would take us 2-3 CPU years on our SGI Indys, and we have not yet attacked it. The most 
time-consuming part of the task is an accurate determination of the transition temperature. 



IV. DETAILS 
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A. Susceptibilities 



The definition ( |1.7| ) of tlie susceptibility is equivalent to 

1 



^ 2N 



(S(0)-S(0))h=o-(S(0))U (4.1) 



or 



X = ^ lim(S(p)* ■ S(p))h=o = lim x(p) , (4.2) 
where is the number of lattice sites and S(p) is the Fourier transform of the spin fields 

S(p)=5:S(x)e^P-^. (4.3) 

X 

To measure the disordered-phase susceptibility in our simulations, we simply use the first 
term of (|]l]). 

For the ordered phase, the subtraction in (|4.1| ) is delicate and we instead use ([1.2|), 
taking long, asymmetric lattices of size L xT xT to get small values of the lowest non-zero 
momentum pmin = 27r/L. In more detail, we find we can get good estimates of X- for fixed 
L by measuring x{p) fo^^ two lowest non-zero momenta in the long direction, 2tt/L and 
4:71 /L, and then extracting X- by fitting to the form 

X-{P)-—T^ 2- (4-4) 

These estimates for X- converge fairly quickly as L is increased. Typical examples of such 
dependence are shown in fig. ^ for x=0.5 and 0.3. For the sake of using all our data, our 



final results for X- in table |T| are the result of a fit of these results for individual L to the 



form^ a + bL ^. There is not too much difference between this and the individual result for 



the largest L. A list of the largest lattice size we use for each x may be found in table IV 



^ From (4^), the difference between x-{Pmin) and the true susceptibihty X- scales as p"^-^^^ ^ L ^ for large 
L. Interpolating x{p) from the two lowest non-zero momenta improves the error to L'^. This is the reason 
for the form of om fit of these interpolated values vs. L. 
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FIG. 12. Measurements of x_ on L xT xT lattices, using the interpolation (O), vs. L. The solid line is 
the best fit of the results to a + bL^'^, and the cross (and smaller error bar) at L—oo shows the extrapolation 
from this fit. Two cases typical of our data are shown: (a) T~6Q, with 63% confidence level; (b) 

x=0.3, r=120, with 18% confidence level. 
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B. Correlation Lengths 



To measure correlation lengths, we first measure the correlation function 

G{r) = (S(0) ■ S(r)) , (4.5) 

where, in practice, we average the right-hand size over all translations and cubic rotations. 
We measure this correlation for the two cases of (1) displacements r that lie along edges of 
the lattice {i.e. in the direction of nearest neighbors), and (2) those that lie along diagonals 
of the lattice. In each case, we also measure the full statistical error matrix of our measure- 



ments, which includes correlated errors between different r. (See sec. [IV E| for our method 
of computing statistical errors.) Then we fit G{r) to the form 



G{r) = a+ 
for the disordered phase and 



-e '''^^^ + (images) 



(4.6) 



G{r) = a. 



-e ^^^^ + (images) 



+ 6- (4.7) 



for the ordered phase. "Images" denotes the similar exponential contributions from nearby 
images of r due to the finite periodicity of the lattice. We first try fitting the above forms to 
all the data points in the relevant direction (edge or diagonal). Then we iteratively throw 
away the smallest r point from our data set until both (1) the confidence level of the fit is 
at least 10%, and (2) r > ^ for all the points being fit. Fig. |1^ shows examples of fits to the 
correlation function. 



Fig. |14| shows the effect of continuing to throw out even more and more short- distance 
points once our criteria are satisfied. Focusing on the fit to ^^*^se^ there is a clear drift of ^ as 
more points are removed, and then there is a plateau that is high compared to our nominal 
value of ^ and its error (the left-most data point ) . One suspects that the value of C, at this 
plateau might be a better estimate than that from our procedure. Unfortunately, we do not 
have a single, universal criterion that would exactly agree with one's best subjective guess 
of ^ for every data set. The drift in values suggests that the systematic errors from fitting 
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may be a bit larger than the (statistical) errors we have assigned to ^. If one strengthened 

the requirement r > ^ for points being fit to r > 1.5^, this drift would change our final 

interpolations for the a; O'^ values of / from ( p.l7| ) to 

ffidgc fdiag 

lim = 1.55(5) , lim ^ = 1.58(6) , (4.8a) 
for extrapolation vs. ^^'^ and 

t-edge fdiag 

£m = 1.60(6) , £m |^ = 1.64(9) , (4.8b) 



for extrapolation vs. This is comparable to the spread of values ( |3.17| ) and consistent 
with our final result (p.. 51). 



The analysis of the correlation lengths for x=0.3 is slightly complicated. Because of the 
large lattice volume required, it takes longer to generate independent configurations than for 
other X. At the same time, we can measure the correlation function G{r) at a larger number 
Nj. = 60 of values of r. If the number of independent configurations is large compared to N^, 
one can estimate the full covariance matrix (see sec. [iV Ej ) for all of these measurements and 



use it to find the correlation length. In our simulations, however, the number of independent 
configurations n generated for x = 0.3 (roughly 70 for the correlation function) is too small 
for this.0 We have chosen to circumvent this issue by simply throwing away many of our r 
values when fitting the correlation function, as we shall describe below. As we shall see, the 
X = 0.3 results do not much affect our final results for ^+/^_; so it is not necessary to make 
a more sophisticated analysis. 

Fig. |15] shows the correlation function along diagonals for x = 0.3 at the transition in the 
ordered phase. Two regions of r plausibly contain the most important information for fitting 



the correlation function: the points at the largest r, which determine 6_ in ( |4.7] ), and the 



points from one to a few correlation lengths, which determine We have chosen to keep 



only the r values marked by diamonds in fig. |I5|: every other point for 6 points at the largest 



^ When n < Nr, for example, the measured covariance matrix will always be singular. 
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x = 0.5, 80x80x80, discordered phase 



x = 0.5, 80x80x80, ordered phase 
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FIG. 13. Examples of fits to the correlation function G{r) in (a) the disordered and (b) ordered phase. 
The scatter of points around the fits is small because of correlations between the points. The dotted vertical 
lines indicate the first point used with the criteria r > ^ (left line, and the value used for the fit) or r > 1.5^ 
(right line). 

r, and every other point for 7 points starting from r ^ ^. We then have roughly five times 
as many independent measurements as r points. We use essentially the same prescription 
for the disordered phase and for the ^jf^*^.[] The results previously quoted in sec. |ITl| were 
obtained using this method. 

We have checked the sensitivity of our results to the choice of which points to keep. If 
we instead simply take 12 to 13 evenly spaced points between r = ^ and Tmax, we obtain 
= 1.42(7) for edges and 1.43(8) for diagonals at x = 0.3, as compared to 1.45(7) 
and 1.41(8) listed in table ^ Due to the large error in our x=0.3 results, this change of 
prescription has negligible effect on the extrapolated ratios of ■C+Z'C-- 



^ Actually, our somewhat arbitrary criteria was to always keep n/Nr as close as possible to 5. As a result, 
in some cases we took only 6 instead of 7 points starting from r « ^. 
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FIG. 14. Examples of the dependence of our fits for ^ on the minimum r kept in the fit. We show ^_ for 
particular runs at /3t for (a) a;=0.6 and (b) x=0.5. ^_ is measured both along edges (squares) and diagonals 
(diamonds). The left-most point of each type is the point chosen for our final value, as described in the 
text. (The left-most point right of the dotted line is chosen if the criteria r > 1.5^ is used instead of r > ^.) 
Confidence levels of the fit are given for each point. 



x=0.3, 120x120x120, ordered phase 
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r 

FIG. 15. Correlation function G{r) along diagonals for x = 0.3 in the ordered phase at the transition. 
The diamonds represent the points actually used for the fitting and are also marked by arrows. 
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C. Finite volume effects 



Table |V| shows the largest lattice sizes we have used in the determination of various 
quantities. For all measurements, every lattice dimension is at least five times the correlation 
length and is typically ten times the correlation length. We shall assess whether finite volume 
errors are significant compared to our statistical errors. Finite size effects should generally 
fall exponentially with increasing lattice size. 



X 




X- 


A 


conversion 


1.0 


403 


160x40^ 


160x40^ 


40 ~ 13.9^?^' 


0.8 




)) 


11 


40 ~ ll.Oe?^' 


0.6 


503 


200x50^ 


11 


50 ~ 9.4^^^'= 


0.5 


803 


240x60^ 


240x60^ 


80 ~ ll.Oe?^' 


0.3 


1203 


240x120^ 


480x120^ 


120 ~ bAC+^' 



TABLE IV. Largest volumes used for various measurements. The last column gives a conversion between 
a typical lattice dimension and the disordered phase correlation length ^^'^sc (^-^yj^jch is larger than ^_). 

In many Monte Carlo applications, the cleanest way to show that finite volume effects 
are negligible for a given lattice volume is to repeat the simulation on smaller and smaller 
volumes until the effects are clearly noticeable. One then extrapolates the finite size cor- 
rections back to the original, large volume. Unfortunately, this procedure is problematical 
in our case. In smaller volumes, the system undergoes transitions between the ordered and 
disordered phases. One can still measure quantities such as C± during a time period between 
transitions, but, as the volume and that transition time decreases, the statistical error of the 
measurement will increase. This degradation of the statistics for smaller volumes obscures 
the finite size effects one would like to measure. 

Instead, we content ourselves with simulating the system for a few different "large" vol- 
umes and checking whether the discrepancies seem consistent with statistical error. Fig. 16 
shows our checks. The data for C+/C_ and look pretty good. The data for x+/X- 
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suggests there might be a systematic effect malting our larger volume measurements slightly 
larger than our lower volume measurements. In any case, we estimate that our finite size 
errors are no larger than our statistical errors. 

We do not have data for x=0.3 but, based on the sizes of the correlation lengths, x=0.3 
on our 120^ lattice should be roughly comparable to x=0.6 on a 30^ lattice. A portion of a 
run on the latter is shown in fig. Not only do we see a transition between the phases, 
but we also see a small spike corresponding to an aborted transition attempt. Whether or 
not such spikes appear in one's data can have a significant effect on the extraction of the 
specific heat. For instance, if the run shown were cut off just before the spike, one would find 
C+ = 2.36(4). If the run were cut-off just after the spike, one would find C+ = 2.51(10) and 
suddenly have drastically revised the error estimate. It is because of this sort of finite-size 
effect that we distinguish our x = 0.3 data as slightly less reliable. 

In a few cases, we have checked the possibility of finite size errors in our determination 
of jSt- (Determining jSt is quite time consuming.) Table |V| shows the results for different 
choices of the transverse size T of the LxT xT lattices we use for determining the transition 
temperature. The values are consistent with each other and, based on the transverse sizes 
of our lattices in units of ^+ (as given by table fVp, we believe that our measurements for 
other X should be reliable as well. 



X 


lattices Pt 


1.0 


L X 20^ 0.157156(12) 
L X 40^ 0.157154(4) 


0.5 


L X 40^ 0.186750(4) 
L X 60^ 0.186751(3) 


0.3 


L X 80^ 0.200362(3) 
L X 120^ 0.2003659(15) 



TABLE V. Dependence of the determination of f3t on transverse lattice size for those x where we 
measured it. 
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FIG. 16. Checks of volume dependence of our ratios for all of our values of x. C'± and S,± were measured 
on an lattice. x+ was measured on a lattice, and x- was extrapolated on L x lattices. Where 
there are two values of T specified, the first is for x- ^-nd the second for x+- For the ^ ratio, the squares and 
diamonds are the results along edges and diagonals respectively. Error bars do not include the uncertainty 
in /Jf 
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FIG. 17. The internal energy vs. Monte Carlo time for an a;=0.6 run on a somewhat smallish lattice 
(30'^). This section of the run shows a transition between the phases and, preceding it, a spike corresponding 
to an unsuccessful transition attempt. 

D. Transition temperature uncertainties 



1. Determining Pt 



As discussed in the introduction, our procedure for determining the transition temper- 
ature is to make multiple runs starting from mixed-phase initial conditions on long, asym- 
metric lattices {L xT xT) and to find the (3 for which the system is equally likely to end up 
in the ordered or disordered phase. A simple, biased random walk model of this process is 
presented in appendix |B], which predicts that the probability P of ending up in the ordered 
rather than disordered phase should have (3 dependence of the form: 



1 + tanh 



(4.9) 



where jSt and A/5 are not determined by the model and A/5 may depend on the transverse 
lattice size and on the Monte Carlo algorithm. An example of our data, and the tanh curve 
that best fits it, is shown in fig. 

At the beginning of each simulation, we obtain the initial condition of fig. ^ by initializing 
one half of the lattice with (3=0 initial conditions, one half with (3=oo initial conditions, and 
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0.1567 0.1568 0.1569 0.157 0.1571 0.1572 0.1573 0.1574 



FIG. 18. An example of the probability of ending in the ordered phase vs. (3 when starting with mixed 
initial conditions. The data is for x = 1.0 (the 4-state Potts model) on a 40 x 20^ lattice. The fit (dashed 



line) to (4.9) has a 75% confidence level. (This is typical of our fits to other data, at smaller x or different 
lattice spacings, where the confidence levels of samples we checked ranged randomly from roughly 45% to 
85%. At smaller x, we generally have somewhat fewer measurements and hence somewhat larger statistical 
uncertainty than shown above.) 

then evolving the two halves independently for roughly 1000 sweeps (depending on x) at 
the desired j3. Only then are the two halves allowed to interact and the interface allowed to 
move, and the entire system is then evolved together. 

Appendix P also discusses why the j3t determined by fitting (|4.9| ) is not necessarily 
correct for finite L, and the model predicts the correction to the true f3t scales like We 



therefore fit Pt for a variety of L and have extrapolated to get the final values of table |T]. 
An example is shown in fig. |19[ Data is shown for two different transverse lattice sizes, 
which have different corrections but extrapolate to consistent L ^ oo limits. 

The simple model of Appendix ^ also predicts that the width Aj3 of the tanh curves 
should scale like 1/L. Though this is not directly relevant to our determination of jSt, it is 



worth checking. Fig. ^ shows a fit of this behavior to the data corresponding to fig. |T^. 
The fit is not very good for the L x 40^ data, with 2% confidence level, and this mediocrity 
is typical of our data at other values of x. In contrast, our fits for /5t typically work fairly 
well. We suspect that the failure of our model for Aj3 may be due to the non-local nature 
of the cluster algorithm. 
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FIG. 19. An example ol the finite L dependence of our determination of /3t. The data is for x^l.O on 
(squares) L x 20^ and (diamonds) L x 40^ lattices. (For L=160, the smaller error bar is the diamond point.) 
The lines are the best fits to the form a + bL^^, and the extrapolations (see table 0) are shown at L—oo. 
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2. Effects on other measurements 



Given our estimates of the uncertainty in the transition temperature (see Table PTl]) , we 
now need to determine the resulting uncertainties in our measurements of C, ^, and %. For 
C and X, we do this by measuring 



a^(r'C) = -iV2((e-e)=^), 

5/3X = -iV[(xe)-(x)(e)] 



(4.10) 
(4.11) 



where iV is the lattice volume, e = EjN is the energy density, and e = (e). We make our 
measurements at the central value of /3f Results are given in table |V|. Our determination 
of these derivatives at small x is not particularly good, but it only needs to be good enough 
to estimate our error due to the uncertainty in /?f For the correlation length, we similarly 
compute the /3 derivatives of the correlation function and of our error matrix. The separate 



sizes of statistical errors and /5t uncertainty errors in our final results is given in table VII 



X 




5/3(/5-'C_) 






1 


2.14(17) X 10^ 


-4.0(6) X 10^ 


1.29(6) X 10^ 


-1.39(20)xl0^ 


0.8 


2.2(3)xl0^ 


-6.5(1.2)xl0^ 


2.55(18)xl0^ 


-1.82(17)xl0^ 


0.6 


4.0(5) X 10^ 


-9.4(1. 9)xl05 


8.7(9) X 10^ 


-6.2(5) xlO^ 


0.5 


6.7(8) X 10^ 


-1.5(4) xlO^ 


2.42(19)xl05 


-1.99(22) xlO^ 


0.3 


4.9(1.2)xl0^ 


-2.3(1. 0)xlO^ 


1.30(24) X 10^ 


-1.21(23)xl0^ 



TABLE VI. Summary of /3 derivatives of C± and x± ^'^ the central value of /3t 
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error / stat. 


error 


X 












edge 


diag. 




1 


0.2 


0.2 


0.3 


0.3 


0.8 


0.2 


0.1 


0.2 


0.2 


0.6 


0.2 


0.3 


0.3 


0.3 


0.5 


0.5 


0.7 


0.8 


0.8 


0.3 


0.4 


0.5 


0.4 


0.4 



TABLE VII. Relative size of error due to uncertainty in j3t compared to the statistical error for the 
ratios whose total errors are given in table ||. The total error is the /3t uncertainty added in quadrature 
with the statistical errors. (For x-i "statistical error" includes the extrapolation L — > oo.) 



E. Statistical errors and algorithms 

Statistical errors were determined by first computing the decorrelation time relevant 
for each quantity. (For example, for C, we use the decorrelation time in the energy.) By 
decorrelation time, we mean a measure of the Monte Carlo time over which configurations 
within a given phase become statistically uncorrelated. Our decorrelation time does not 
measure the mixing time between phases, since we run in large enough volumes that there 
are no transitions between phases in our simulations. 

Consider first estimating the error of quantities which are determined by simple ensemble 
averages of some quantity A {e.g. the susceptibilities or the energy density e). We then 
use the integrated decorrelation time defined by [[l^j 

1 +00 

rint = - E C{A-t). (4.12) 

^ t=-oo 

C{A] t) is the auto-correlation function, estimated for a sample of n measurements as 



1 



n—t 



■Y.{A,-A){A,^,-A) 



n — t 

C{A- 1) = 2 ' (4-13) 
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where A is the sample average of A and 



1 " _ 

4 = -E(A,-Af. (4.14) 



The error in A is then estimated as 



ETT(A) = aAJ^ (4.15) 
V n 



In practice, the sum in the definition ([4.1 2|) must be cut off, because the statistical error 



in C(A;t) itself becomes large for large t. Our criteria is to cut off the sum when C{A;t) 
drops below 0.05. The value 0.05 was chosen to give reasonable agreement with the binning 
method described next. Reducing 0.05 to 0.01 would not appreciably change our results. 

Alternatively, one can estimate errors by binning the sequence of measurements Ai {i = 
1, n) into bins of some size Tbin and averaging the data over each bin to obtain a new sequence 
Ak {k = 1, Abin), where A^bm = f^/'Thin- The error is then estimated by the variance of this 
new sequence: 

[Err(:4)]2 = ^ UA') - {Af] . (4.16) 

UVbin — i ) 

Fig. pi] shows a typical example of the result vs. bin size. The error stabilizes at large bin 



sizes, as it should, and roughly agrees with the integrated decorrelation time method, which 
is also shown in the figure. 

To calculate correlation lengths, we need the full correlated error matrix (the covariance 
matrix) of the correlation functions G{r). We compute this using the binning method. 
Binning (with fixed bin size) has the advantage that the resulting covariance matrix is 
positive definite. The covariance matrix is given by 

= j^-^A{Q{n)Q{rj)) - {Q{ri)){Q{r,))] , (4.17) 

l^Vbin — ij 

where ^(r) = S(r) ■ S(0) averaged over translations and cubic rotations. In order to have a 
simple universal criteria for what size of bin to use, we have surveyed a variety of examples 
and found that a bin size of lOrj works well, where Td is the maximum (over the range of r 
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FIG. 21. Example of statistical error estimates vs. bin size. The data is from measurements of x+ at 
X ~ 0.5. The vertical axis has been arbitrarily normalized to be 1.0 for Thin — 2Trf. The solid horizontal line 
is the independent error estimate using the integrated decorrelation time. The dotted line is what the same 



estimate would be if we had cut off the sum ( 4.12 ) when C{A; t) dropped below 0.01 rather than 0.05. 



used for the fit) of tlie integrated decorrelation time for the G{r). Fig. p^ a shows a typical 
example of the dependence of the statistical error of correlation length on bin size, and our 
particular choice of bin size can be seen to be adequately large. 

The error estimate for the correlation lengths are determined by a standard chi-square 
analysis. 

One case that requires a slightly different approach is the specific heat, which is measured 
from the variance of the energy density. In this case, we bin the energy density as above 
but don't average it over each bin. Then we apply the jackknife procedure on the bins: 
we compute the specific heat after throwing away the Tbin energy density measurements 
corresponding to one of the bins. The A^bin possibilities for which bin to have thrown out 
gives us A'bin separate measurements of the specific heat. The variance of these measurements 
is used for the error estimate on the specific heat. By surveying the dependence on bin size 
for all are data, we have found that a good choice is Tbin = lOrd where is the integrated 
decorrelation time for the energy. An example of the dependence on bin size is shown in 
Fig. Hb. 
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FIG. 22. Examples of the dependence of statistical error on bin size for (a) and (b) The 

vertical axis has been arbitrarily normalized to be 1.0 for r^in = ^Td- The diamond marks the actual bin 
size used according to our criteria. 

A similar approach is used for X- and all (3 derivatives [e.g. dpx±)- 

We have checked that running simulations with a pure heatbath algorithm, or with a 

pure cluster update algorithm, give results that are statistically consistent with each other 

and with our interlacing of the two algorithms.^ 

This work was supported by the U.S. Department of Energy, grants DE-FG06-91ER40614 
and DE-FG03-96ER40956. We thank Joseph Rudnick, Michael Fisher, Joan Adler, and 
David Wright for useful conversations, Peter Ungar for explaining the solution to the Gam- 
bler's Ruin problem. Marcel Den Nijs for pointing out that the cross-over exponents are 



^ In fact, such checks led us to discover for ourselves the problems of naively using simple random- number 
generators. We had initially used a 32-bit congruence algorithm with a period of . This period is too 
small for the length of some of our simulations, manifesting in inconsistency between the heatbath and 
cluster algorithms and producing correlations which could clearly be seen in the very-long-time tail of the 
energy auto-correlation function. We had to switch to a generator with period 2^*^. 
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determined by well-known Ising model critical exponents, and, most of all, Steve Sharpe 
and Larry Yaffe, who gave too many valuable suggestions to enumerate. 



APPENDIX A: CLUSTER UPDATE ALGORITHM 

The two of us have used different versions of a cluster update algorithm, one which grows 
and flips a single cluster at a time, and another which grows simultaneous clusters across 
the entire lattice. In any case, they are simple generalizations of the algorithms ||T8| , |l3| used 
in other spin systems. Write the Hamiltonian as 

(3H = Y.hiS,,S,) (Al) 

where h{Si, Sj) is the nearest-neighbor interaction. 

One step of the first version of the algorithm can be summarized as follows, (a) Randomly 
choose one of the five order-2 elements R {R^ = 1) of the internal symmetry group D^. (b) 
Randomly choose a lattice site x as the first point to include in the cluster c, and mark the 
site, (c) One at a time, visit all new links (xy) connecting a; G c to its nearest neighbors y. 
For each link visited, check if y is already in c, and if not then adjoin it to c with probability 
P(Si., Sy), where 

P(S,, S,) = 1 - exp {min[0, h{S,, S,) - h{RS,, Sy)]} . (A2) 

A newly included y should be marked, (d) Repeat step (c) until no new sites are added to 
the cluster, (e) Flip all the spins in the resulting cluster: S ^ RS. 

Randomness of the choices in steps (a) and (b) is inessential: one only needs to vary the 
choices enough to give reasonably efficient ergodicity. If (a) is restricted to a random choice 
between just the two symmetries corresponding to s —s and t ^ —t respectively, then 
this algorithm is equivalent to the Ashkin- Teller cluster algorithm described in ref. [|1^]. (See 



also ref. ^|. 



One step of the second version of the algorithm is as follows, (a) Randomly choose an 
R as above, (b) Visit every link {xy) in the lattice and mark it with probability P(Sx>, S^^) 
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given by ( |A2|) . (c) Identify all the disconnected clusters of sites connected by marked links, 
and flip each such cluster with probability 1/2. 



APPENDIX B: THE GAMBLER'S RUIN PROBLEM 

In section ||, we discussed how we determine the transition temperature by splitting an 
asymmetric L x T x T lattice in half along the longitudinal direction L and starting the 
two halves in the ordered and disordered phases respectively. In this appendix, we explain 
a simple model for the probability of the system evolving into one phase over the other. 

At any time, let Zi and Z2 be the (transversely averaged) locations of the two domain 
walls, and let z be the separation between them as measured through the disordered phase 
(pi L — z as measured through the ordered phase). When /3=/3t, z should random walk 
in Monte Carlo time. When (3 deviates slightly from /5t, there will be a slight bias in this 
random walk proportional to (3—(3t. If we model this biased random walk as taking fixed 
steps in z with probability 

prob(2 ^ 2+1) = p = 1(1 + e) , (Bl) 
Y>ioh{z ^ z-l) = q = \{l - t) , (B2) 

where e oc (3—(3t is small, then we have a special case of the Gambler's Ruin problem. The 
Gambler's Ruin problem is that you start with z dollars in your pocket and you play some 
casino game over and over again until you either go broke or accumulate your goal of L 
dollars. What is the probability you can afford the taxi home? 

To solve it, let P{z) be the probability of winning if you start at z. Then, by considering 
one step, one obtains the difference equation 

P{z)=pP{z + l) + qP{z-l) , (B3) 

and the boundary conditions are 

P(0) = 0, P(L) = 1. (B4) 
41 



The solution is 




(B5) 



In our case (small e and large L), this becomes 



1 -e 



(B6) 



1-e 



-2eL ■ 



Putting in our initial condition z = L/2, we find that the curve in fig. ^ is a tanh curve: 



P{L/2) - i[l + tanh(ieL)]. 



(B7) 



We don't know a priori the proportionality constant between our parameter e and f3 — (3t, 
but we can parametrize the curve as 



and determine jSt and A/3 by fitting to our Monte Carlo results. We haven't attempted to 
model the dependence of e on transverse size T. However, for fixed transverse size the model 
([B7D predicts 



As mentioned in sec. 0, there are systematic errors which make the center of the tanh 
curves shift as L is increased. One example would occur if the two halves aren't adequately 
equilibrated at the inverse temperature f3 before allowing the domain walls to evolve. Then, 
even exactly at (3t, there might be some systematic bias to the initial motion of the walls 
until thermal equilibrium is reached. This could be modeled by simply starting z with some 
systematic initial offset, i.e. z = L/2 + a. Secondly, even if proper equilibrium is reached, 
the "stickiness" of the domain walls need not be the same in the two phases. For example, 
we have only been discussing the average longitudinal coordinates zi and Z2 of the domain 
walls. In fact, the domain walls have transverse excitations and might have slightly longer, 
thinner fingers reaching out into one phase than into the other. The separation at which 




(B8) 



Ap oc 1/L. 



(B9) 
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they first touch each other might be larger in one phase than in the other. Such an effect 
could be modeled by considering the effective end-points of the game to be slightly and 
asymmetrically different from z=0 and z=L, now being z=bi and z=L — 62- By a shift 
of coordinate, this can again be considered as the problem of starting slightly away from 
z = L/2. 

So a model for the systematic error is to ask what happens if the probability is really 
P{L/2 + a) from (|B^ ), for some a, but we nonetheless tried to extract a value of /3t fitting 
the form P{L/2) for (|B7| ) to the result. One easily finds that the systematic error in /3t 
then scales as for large L (assuming fixed transverse dimension and hence a fixed 

proportionality between e and P — Pt)- 



APPENDIX C: BCC LATTICES 



We have made a few simulations on BCC lattices, but not enough to extract any a;— >0 
limits. Our results are presented in table |V114 The measurements of C± were made on 40^ 
lattices; jSt was measured on lattices as big as 160 x 20^ for x=1.0 and 80 x 40^ for x=0.6. 
When we refer to an Li x L2 x L3 BCC lattice, we mean one with L1L2L3 unit cells and so 
= 2L1L2L3 sites. All of our lattices are helical. Defining helical boundary conditions on 
a BCC lattice is perhaps non-standard, and we explain it below. 



X (3t 




c+ 






e_ 






BCC: 








1 0.113752(5) 


-2.344(4) 


1.34(24) 


0.112(3) 




-5.808(5) 


11.9(3) 




0.6 0.130291(6) 


-2.1975(8) 


1.88(4) 


0.104(4) 




-4.392(7) 


18.1(5) 
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TABLE VIII. Same as table III but for BCC lattices. Densities are given in units of per lattice point, 
which for BCC differs from per unit cell. 



Helical boundary conditions for BCC lattices 

To motivate our definition of lielical boundary conditions for BCC lattices, we first 
briefly review the situation for simple cubic lattices and, for simplicity of presentation. 



will first consider two-dimensional examples. Fig. ^ shows an infinite simple cubic lattice 
in two dimensions. A finite-volume lattice with periodic boundary conditions corresponds 
to restricting the system to the sites inside a box, such as shown in fig. [23|a, and then 
identifying opposite edges of the box. The same volume with helical boundary conditions 



corresponds instead to the box in fig. P3|b, again with opposite edges identified. Helical 
boundary conditions retain translation invariance. Their advantage is that the sites can 
be numbered, as shown in fig. |23|b, in such a way that the offset between one site and its 
neighbor in a given direction is always a fixed number modulo N, independent of the site 
chosen, where N is the total number of sites. For the two-dimensional lattice shown, the 
offsets An corresponding to the four directions are 

An = ±1, ±Li mod N=LiL2 . (CI) 

In the three-dimensional case, they are 

An = ±1, ±Li, ±LiL2 mod N=LiL2L^ . (C2) 

The nature of these offsets is an advantage because, if the lattice is represented as a linear 
array in memory, it makes indexing neighbors of sites quicker and easier than for periodic 
lattices. (|C^ ) is the definition used in this paper for an Li x L2 x L3 helical simple cubic 
lattice. 

Fig. 0a shows a somewhat similar box drawn for a BCC lattice.0 Unfortunately, there 



^ In two dimensions, an infinite BCC lattice is of course equivalent to a simple cubic lattice. We are 
discussing it just as a visual aid for generalizing to three-dimensional BCC lattices. 
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FIG. 23. Ways to choose a finite volume box from the infinite plane and impose periodicity correspond- 
ing to (a) normal periodic boundary conditions, and (b) helical boundary conditions. 
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FIG. 24. Ways to try to choose a finite volume box for a BCC lattice, (a) lacks the advantages of the 
helical simple cubic case; (b) is our definition of a helical BCC lattice. 
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is no way to number the sites in the box so that the offset An between neighbors is the 
same for all sites. However, if one instead draws the box as in fig. p4|b, such numberings are 
possible. For two- dimensions, the numbering of fig. corresponds to 

An = ±Li, ±(Li - 1) mod N=2LiL2 , (C3) 

For three dimensions, where there are 8 nearest-neighbor directions, the generalization is 

An = iLiLa, ±(^1^2 - 1), ±(^1^2 - Li), ±(^1^2 - Li - 1) mod N=2LiL2L3 . (C4) 

This is our definition of an Li x L2 x L3 helical BCC lattice. 

APPENDIX D: A MODEL OF RATIOS AND CROSS-OVER EXPONENTS FOR 

SMALL a 

1. Overall scaling 

In this appendix, we elaborate on the assertions of sec. |111 A| about how small x needs 



to be if the Ising critical exponent a is formally considered small. We begin by reviewing 
the case of non-small a in slightly more detail. For the sake of specificity, we focus on the 
behavior of the correlation lengths The correlation functions G±{R) = {s{R)s{0))± will 
have the following form of universal cross-over scaling function near the transition for small 



x: 



G±{R) = h^yG^ipyH, W^x, b-^R) , (Dl) 

where b is the arbitrary renormalization distance scale and 

y = d-2 + r], yt = l/iy, y^ = a/iy (D2) 

are the anomalous dimensions of G±, t, and x. We have ignored the effects of irrelevant 
operators, t is a scaling field corresponding to the reduced temperature. In this language. 
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the derivation of the relationship between ^± and x can be made by first choosing b = x ^/^"^ 
and writing 

G±{R) = xy/y^G±{x-y'/yH, i, x^^y^R) . (ds) 

As we vary t, the transition must occur for some definite value Tq of the first-argument on 
the left-hand side: 

G±{R) = xy/y^G±{To, 1, x^/y^R) . (D4) 



By dimensional analysis, the long-distance behavior of the right-hand side of must have 
the form 

G±(ro,l,r)~e"'^/^±, (D5) 

and so 

e± ~ A±x-^/y=^ , (D6a) 

^- A_ ^ ' 

Now we can discuss the case where a is small or zero. ( pi|) was a special case of a more 
general situation where x doesn't scale as a power law. To be more general, replace 

by^x x{b) , (D7) 

where x{b) is the solution to some renormalization-group equation. The non-small a case 
corresponded to 

bdbX = y^x + 0{x^) = -x + 0{x^), (D8) 

The form ( pi|) corresponds to ignoring the O(x^) correction for small x, giving a simple, 
power-law solution for x. When a = 0, on the other hand, the O(x^) terms in ( pi|) become 
essential: 

bdbX = cx^ + 0{x^) , (D9) 
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where we assume c > 0. (c < would lead to a second-order Ising phase transition for small 
X > and so is excluded if we assume the phase transition remains first-order.) Ignoring 
the 0{x^) terms, this corresponds to our original scaling function (pi|) with 



b^^x \ , . (DIO) 

X ^ — c In 

For very small x, we now want to choose 

b ~ e^/"^ , (Dll) 



which gives ^± ~ e^^^^ as discussed in sec. |1II A| . However, just as in ( P6b|) , one will still get 



Now we're ready to model the case where a is arbitrarily small but non-zero. Clearly we 
do not want to ignore the O(x^) term in the renormalization group equation, so we take 

bdbX = - x + cx'^ + 0{x^). (D13) 

Ignoring higher-order terms, the solution is 

which gives 

\ vc ' 

which interpolates between the previous cases and has the properties summarized in 



sec. 



Ill A| . But, just as in the previous cases, again 

|l: = 4i^a;0_ pig) 

So the ratio is insensitive to the crossover between exponential and power-law dependence 
of ^± on x~^ . 
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2. Corrections to scaling 



To address corrections to scaling, supplement ( |D1| ) by the most important irrelevant 
operator, whose coefficient we shall call 

G±{R) = h-yG±{WH, W^x, b-^R, b-'^z) . (D17) 

Now choosing b as before gives 

G±{R) = xy/y^G±{x-y'/yH, i, x^^^^r, x'^iy^z) . (dis) 

Quantities like A± now depend on x'^l'^'^z but, for small can be Taylor expanded: 

|i = ^ = 0(x°) [l + 0(x"/^-)] . (D19) 

For small or zero a, we should make the replacement (P7|) as before, and turn to the 
second-order RG equation (|D13|) .P^ One easily finds that correction-to-scaling laws such as 
( p.l6| ), which do not explicitly depend on a, remain valid. 



Some readers may be more familiar with thinking of the RG flow in the cubic anisotropy model , 
which is in the same universality class, and in terms of the e expansion. The potential energy in that model is 
of the form + +(/)2). Roughly speaking, z here corresponds to v and x to a linear combination 

of —u and v. 
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There will in general also be a term dxz on the right-hand size of (D13) but, because z is irrelevant, 



this term quickly becomes negligible as h is increased and does not affect any of our conclusions. 
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